American mink (Neovison vison) space use in Illinois (data from Ahlers et al. 2015)

packages<-c("adehabitatHR","data.table","ggfortify","grid","move","moveVis","OpenStreetMap","pbapply","plotly","rgdal","sp","tidyverse","viridis")
sapply(packages, require, character.only=T)

 adehabitatHR    data.table     ggfortify          grid          move       moveVis OpenStreetMap       pbapply 
         TRUE          TRUE          TRUE          TRUE          TRUE          TRUE          TRUE          TRUE 
       plotly         rgdal            sp     tidyverse       viridis 
         TRUE          TRUE          TRUE          TRUE          TRUE 

American Mink Dataset

This dataset looked at the American Mink home ranges in Illinois where there habitat has been highly altered.They studied 20 individuals for 6 years (2007-2012) and recorded present. The years in which they studied the minks they had record-breaking drought and flood conditions. They used mark-recapture methods and telemetry to track the individuals.

For this assignment I looked at individuals 1 , 10 and 11

Individual 18: 192 locations

Individual 19: 113 locations

Individual 7: 117 locations

American_Mink <- read_csv("American_Mink.csv")

qaqc_plot <- ggplot() + geom_point(data=American_Mink, 
                                   aes(utm_easting,utm_northing,
                                       color=individual_local_identifier)) +
                        labs(x="Easting", y="Northing") +
                        guides(color=guide_legend("Identifier"))

ggplotly(qaqc_plot)

lapply(split(American_Mink, American_Mink$individual_local_identifier), 
       function(x)write.csv(x, file = paste(x$individual_local_identifier[1],".csv"), row.names = FALSE))

$`7`
NULL

$`18`
NULL

$`19`
NULL

files <- c("18 .csv","19 .csv", "7 .csv")

utm_points <- cbind(American_Mink$utm_easting, American_Mink$utm_northing)
utm_locations <- SpatialPoints(utm_points, 
                 proj4string=CRS("+proj=utm +zone=16 +datum=WGS84"))
proj_lat.lon <- as.data.frame(spTransform(
                utm_locations, CRS("+proj=longlat +datum=WGS84")))
colnames(proj_lat.lon) <- c("x","y")
raster <- openmap(c(max(proj_lat.lon$y)+0.01, min(proj_lat.lon$x)-0.01), 
                  c(min(proj_lat.lon$y)-0.01, max(proj_lat.lon$x)+0.01), 
                  type = "bing")
raster_utm <- openproj(raster, 
              projection = "+proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs")

autoplot(raster_utm, expand = TRUE) + theme_bw() +
  theme(legend.position="bottom") +
  theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
  geom_point(data=American_Mink, aes(utm_easting,utm_northing,
             color=factor(individual_local_identifier)), size = 3, alpha = 0.8) +
  theme(axis.title = element_text(face="bold")) + labs(x="Easting",
        y="Northing") + guides(color=guide_legend("Identifier"))

Minimum Convex Polygon

library(pbapply)
mcp_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm_easting)
  y <- as.data.frame(data$utm_northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data , proj4string = CRS("+proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  mcp.out <- mcp(xy, percent=100, unout="ha")
  mcp.points <- cbind((data.frame(xy)),data$individual_local_identifier)
  colnames(mcp.points) <- c("x","y", "identifier")
  mcp.poly <- fortify(mcp.out, region = "id")
  units <- grid.text(paste(round(mcp.out@data$area,2),"ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  mcp.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=mcp.poly, aes(x=mcp.poly$long, y=mcp.poly$lat), alpha=0.8) +
    geom_point(data=mcp.points, aes(x=x, y=y)) + 
    labs(x="Easting (m)", y="Northing (m)", title=mcp.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  mcp.plot
}

pblapply(files, mcp_raster)

  |                                                                                  | 0 % ~calculating  
  |============================                                                      | 33% ~01s          
  |=======================================================                           | 67% ~00s          
  |==================================================================================| 100% elapsed=01s  
[[1]]

[[2]]

[[3]]

Kernel-Density Estimation

kde_raster <- function(filename){
  data <- read.csv(file = filename)
  x <- as.data.frame(data$utm_easting)
  y <- as.data.frame(data$utm_northing)
  xy <- c(x,y)
  data.proj <- SpatialPointsDataFrame(xy,data, proj4string = CRS("+proj=utm +zone=16 +ellps=WGS84 +units=m +no_defs"))
  xy <- SpatialPoints(data.proj@coords)
  kde<-kernelUD(xy, h="href", kern="bivnorm", grid=100)
  ver <- getverticeshr(kde, 75)
  kde.points <- cbind((data.frame(data.proj@coords)),data$individual_local_identifier)
  colnames(kde.points) <- c("x","y","identifier")
  kde.poly <- fortify(ver, region = "id")
  units <- grid.text(paste(round(ver$area,2)," ha"), x=0.85,  y=0.95,
                     gp=gpar(fontface=4, col="white", cex=0.9), draw = FALSE)
  kde.plot <- autoplot(raster_utm, expand = TRUE) + theme_bw() + theme(legend.position="none") +
    theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) +
    geom_polygon(data=kde.poly, aes(x=kde.poly$long, y=kde.poly$lat), alpha = 0.8) +
    geom_point(data=kde.points, aes(x=x, y=y)) +
    labs(x="Easting (m)", y="Northing (m)", title=kde.points$identifier) +
    theme(legend.position="none", plot.title = element_text(face = "bold", hjust = 0.5)) + 
    annotation_custom(units)
  kde.plot
}

pblapply(files, kde_raster)

  |                                                                                  | 0 % ~calculating  
  |============================                                                      | 33% ~00s          
  |=======================================================                           | 67% ~00s          
  |==================================================================================| 100% elapsed=00s  
[[1]]

[[2]]

[[3]]